In [1]:
# Importing
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
%load_ext autoreload
%autoreload 2
import pandas as pd
import numpy as np
import os
import islets
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'
import islets._manuscript_functions as mf
import statsmodels.api as sm
In [2]:
pathToFile = "/data/Sandra/2022/3xTg_AD_w_legs.csv"
print (open(pathToFile).read())
In [3]:
data = pd.read_csv(pathToFile)
data["legs"] = data["legs"].apply(eval)
data
Out[3]:
In [4]:
legList = list(np.unique(sum(data["legs"].apply(list),[])))
In [5]:
legList
Out[5]:
In [6]:
legColors = dict(zip(legList, plt.cm.Dark2.colors))
# lets also plot them to see if python understands what we wrote and if we like them
_, ax = plt.subplots(1,1,figsize = (4,len(legColors)/3))
for il,leg in enumerate(legColors):
ax.axhline(il,color=legColors[leg])
ax.text(0,il,leg+" ", va="center", ha="right", color=legColors[leg])
mf.mystyle_axes(ax)
ax.set_title("leg colors")
Out[6]:
In [7]:
# variable of interest#
voi = "genotype"
In [8]:
data = data.sort_values([voi])
In [9]:
data[voi].unique()
Out[9]:
In [10]:
# Here you specify the order in which they appear
varList = ['wt', 'tg']
In [11]:
varColorCode = dict(zip(varList, ["black","Grey"]))
# lets also plot them to see if python understands what we wrote and if we like them
_, ax = plt.subplots(1,1,figsize = (4,len(varColorCode)/3))
for il,var in enumerate(varColorCode):
ax.axhline(il,color=varColorCode[var],lw=2)
ax.text(0,il,var+" ", va="center", ha="right", color=varColorCode[var])
mf.mystyle_axes(ax)
Data preparation¶
In [12]:
nc = 8
nr = int(np.ceil(len(data)/nc))
fig, axs = plt.subplots(nr,nc,figsize=(3*nc,3*nr))
ia = 0
for i in data.index:
path = data.loc[i,"path"]
## here we load the regions if we didn't ye
if "regions" in data.columns and isinstance(data.loc[i,"regions"], islets.Regions):
regions = data.loc[i,"regions"]
else:
regions = islets.load_regions(path)
data.loc[i,"regions"] = [regions]
## here we plot the regions
ax = axs.flat[ia]
regions.plotEdges(ax=ax, scaleFontSize=8)
ia += 1
ax.set_xticks([])
ax.set_yticks([])
var = data.loc[i,voi]
ax.text(0,-.01,f"{var}",size=8,va="top",transform=ax.transAxes, color=varColorCode[var])
ax.text(1,-.01,f"exp: {i}",size=8,va="top",transform=ax.transAxes, color=varColorCode[var],ha="right")
# remove unused axes
for ax in axs.flat[ia:]:
ax.remove()
In [13]:
for i in data.index:
# construct events_path from the roi path
path = data.loc[i,"path"]
events_path = path.replace("_rois.pkl","_auto_events.csv")
# regions are in our table now
regions = data.loc[i,"regions"]
# reading events each time from the file is very fast
Events = pd.read_csv(events_path, index_col=0)
# plot events
fig, axs = mf.plot_events(Events, regions)
ax = axs[0]
var = data.loc[i,voi]
legs = data.loc[i,"legs"]
ax.text(0,1.3,f"{voi}: {var} exp: {i}",size=8,va="top",transform=ax.transAxes)
for leg in legs:
mf.emphasize_region(ax,legs[leg],ax.get_ylim(), color=legColors[leg])
ax.text(legs[leg][0],ax.get_ylim()[1]*.9," "+leg,va="top",color=legColors[leg])
ax.set_ylim(ax.get_ylim()[0]*.8,None)
plt.show()
In [14]:
hwBinEdges = np.geomspace(start = .1, stop = 1000, num = 100)
hwBinCenters = np.sqrt( hwBinEdges[1:] * hwBinEdges[:-1])
fig, axs = plt.subplots(1,2,figsize=(14,6))
axs[0].set_title("pdf")
axs[1].set_title("cdf")
for ax in axs:
ax.set_xscale("log")
for i in data.index:
# construct events_path from the roi path
path = data.loc[i,"path"]
events_path = path.replace("_rois.pkl","_auto_events.csv")
# reading events each time from the file is very fast
Events = pd.read_csv(events_path, index_col=0)
x = Events["halfwidth"].values.copy()
h = np.histogram(x, hwBinEdges)[0]
h = h/h.sum()
axs[0].plot(hwBinCenters, h, label=f"exp: {i}")
x.sort()
axs[1].plot(x, np.linspace(0,1,len(x)), label=f"exp: {i}")
axs[0].set_ylim(0,)
axs[0].set_yticks([])
axs[0].set_ylabel("arbitrary units")
axs[1].set_ylim(0,1)
for ax in axs:
ax.set_xlim(hwBinCenters[[0,-1]])
axs[1].legend(loc=4)
Out[14]:
Analysis¶
In [15]:
hwRegions = {
# "ultrashort": (0.2,1.5),
"short": (0.5,20),
# "long": (20,200),
}
In [16]:
# we again iterate over rows of data, this time collecting the important objects for each experiment:
allEpmpr = []
minEvents = 10
# ...
for i in data.index:
# ... import events, ...
path = data.loc[i,"path"]
events_path = path.replace("_rois.pkl","_auto_events.csv")
Events = pd.read_csv(events_path, index_col=0)
# define legs
islets.EventDistillery.define_legs(Events, data.loc[i,"legs"])
# ... calculate events per min per active roi, ...
evpmpar = mf.get_events_per_min_per_nrois( # <-- see Step3 for an explanation what this does
Events,
hwRegions,
minEvents=minEvents,
reduceRois=2,
Nbootstrap=30
)
# ... we label result with the experiment index
expLabel = "exp"+str(i)
evpmpar["experiment"] = expLabel
evpmpar[voi] = data.loc[i,voi]
## and collect it into the big list
allEpmpr += [evpmpar]
In [17]:
allEpmpr = pd.concat(allEpmpr,ignore_index=True)
In [18]:
# check one by eye to make sure it makes sense
allEpmpr.query("experiment=='exp0'").sort_values(["genotype","leg"])
Out[18]:
In [19]:
# allEpmpr.to_csv("events_per_min_per_active_roi.csv")
In [20]:
#### Details
#mf.check_roi_heterogeneity(Events.query("halfwidth<20"),regions, boxplot_kwargs={"showfliers":False});
In [21]:
exps = list(allEpmpr["experiment"].unique())
In [22]:
jitter = .01 # how much you wish to stagger the points in the plot below
allEpmpr["plot_pos"] = allEpmpr["leg"].apply(lambda xi: legList.index(xi)) + allEpmpr["experiment"].apply(lambda xi: exps.index(xi)-(len(exps)-1)/2)*jitter
In [23]:
rescale = 100 # rescale to values you can change it if you wish
fig, axs = plt.subplots(len(varList),len(hwRegions),figsize=(1+len(hwRegions)*len(legColors)*.5,3*len(varList)), squeeze=False, sharex="col", sharey="col")
# plt.subplots_adjust(wspace=.5, top = 1-.2/axs.shape[1], bottom = .3/axs.shape[1])
if rescale>1:
fig.suptitle(f"events per minute per {rescale} active rois*")
else:
fig.suptitle(f"events per minute per active roi*\n")
legNames = ["6","S","S+5","S+10","S+20","S-2"]
for i_col,hwr in enumerate(hwRegions):
for i_row,var in enumerate(varList):
ax = axs[i_row, i_col]
for exp,df in allEpmpr.query(f"hw_region=='{hwr}' and {voi}=='{var}'").groupby("experiment"):
ax.plot(df['plot_pos'],rescale*df['epmpr'], color="grey",lw=.7)
for leg in legs:
df = allEpmpr.query(f"hw_region=='{hwr}' and {voi}=='{var}' and leg=='{leg}'").copy()
ax.errorbar(df['plot_pos'],df["epmpr"]*rescale,df["epmpr_std"]*rescale,color=legColors[leg], marker="o",ls="none", mec="k")
if i_row==0:
ax.text(0,1.03,hwr, transform=ax.transAxes)
ax.text(1,.95,f"{voi}: {var} ", transform=ax.transAxes,ha="right")
ax.set_xlim(-.5,int(allEpmpr["plot_pos"].max())+.5)
# ax.set_ylim(0,)
# set the ticks:
for ax in axs.flat:
xmax = ax.get_xlim()[1]
ax.set_xticks(np.arange(int(xmax)+1))
ax.set_xticklabels(legNames, rotation=0)
ax = axs[0,0]
ax.text(1,1.03,r"$\tau_{1/2}<%g$"%hwRegions["short"][1],ha="right",transform=ax.transAxes)
fig.text(0,-.01,f"*active rois are rois that have at least {minEvents} events detected",va="top", )
Out[23]:
In [24]:
from tqdm.notebook import tqdm
In [25]:
# we again iterate over rows of data, this time collecting the important objects for each experiment:
allHWs = []
minEvents = 10
for i in tqdm(data.index):
# ... import events, ...
path = data.loc[i,"path"]
events_path = path.replace("_rois.pkl","_auto_events.csv")
Events = pd.read_csv(events_path, index_col=0)
# define legs
legs = data.loc[i,"legs"]
islets.EventDistillery.define_legs(Events, legs)
Events = Events.dropna().copy()
for hwr in hwRegions:
hw0,hw1 = hwRegions[hwr]
for leg in legList:
evs = Events.query(f"halfwidth<{hw1} and halfwidth>={hw0} and leg=='{leg}'")
x = [dfr["halfwidth"].median() for roi,dfr in evs.groupby("roi") if len(dfr)>minEvents]
allHWs += [{
"experiment": "exp"+str(i),
"leg": leg,
voi: data.loc[i,voi],
"halfwidth": np.mean(x),
"halfwidth sem": np.std(x)/np.sqrt(len(x)-1),
"hw_region": hwr,
}]
In [26]:
allHWs = pd.DataFrame(allHWs)
In [27]:
allHWs.query("experiment=='exp0'")
Out[27]:
In [28]:
jitter = .01 # how much you wish to stagger the points in the plot below
allHWs["plot_pos"] = allHWs["leg"].apply(lambda xi: legList.index(xi)) + allHWs["experiment"].apply(lambda xi: exps.index(xi)-(len(exps)-1)/2)*jitter
In [29]:
fig, axs = plt.subplots(len(varList),len(hwRegions),figsize=(1+len(hwRegions)*len(legColors)*.5,3*len(varList)), squeeze=False, sharex="col", sharey="col")
# plt.subplots_adjust(wspace=.5, top = 1-.2/axs.shape[1], bottom = .3/axs.shape[1])
fig.suptitle(f"median active rois*")
legNames = ["6","S","S+5","S+10","S+20","S-2"]
for i_col,hwr in enumerate(hwRegions):
for i_row,var in enumerate(varList):
ax = axs[i_row, i_col]
for exp,df in allHWs.query(f"hw_region=='{hwr}' and {voi}=='{var}'").groupby("experiment"):
df = df.copy().dropna()
ax.plot(df['plot_pos'],df['halfwidth'], color="grey",lw=.7)
for leg in legs:
df = allHWs.query(f"hw_region=='{hwr}' and {voi}=='{var}' and leg=='{leg}'").copy()
ax.errorbar(df['plot_pos'], df["halfwidth"], df["halfwidth sem"], color=legColors[leg], marker="o",ls="none", mec="k")
if i_row==0:
ax.text(0,1.03,hwr, transform=ax.transAxes)
ax.text(1,.95,f"{voi}: {var} ", transform=ax.transAxes,ha="right")
ax.set_xlim(-.5,int(allHWs["plot_pos"].max())+.5)
# ax.set_ylim(0,)
# set the ticks:
for ax in axs.flat:
xmax = ax.get_xlim()[1]
ax.set_xticks(np.arange(int(xmax)+1))
ax.set_xticklabels(legNames, rotation=0)
ax = axs[0,0]
ax.text(1,1.03,r"$\tau_{1/2}<%g$"%hwRegions["short"][1],ha="right",transform=ax.transAxes)
fig.text(0,-.01,f"*active rois are rois that have at least {minEvents} events detected",va="top", )
Out[29]:
